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Abstract. An iterative scheme based on the kernel polynomial method is devised 
for the efficient computation of the one-body density matrix of weakly interacting 
Bose gases within Bogoliubov theory. This scheme is used to analyze the coherence 
properties of disordered bosons in one and two dimensions. In the one-dimensional 
geometry, we examine the quantum phase transition between superfluid and Bose 
glass at weak interactions, and we recover the scaling of the phase boundary that was 
characterized using a direct spectral approach by Fontanesi et al. [Phys. Rev. A 81, 
053603 (2010)]. The kernel polynomial scheme is also used to study the disorder- 
induced condensate depletion in the two-dimensional geometry. Our approach paves 
the way for an analysis of coherence properties of Bose gases across the superfluid- 
insulator transition in two and three dimensions. 



1. Introduction 

More than twenty years ago, the discovery that superfluidity may be suppressed in ^He 
adsorbed on porous media [1, 2, 3, 4] triggered investigations into the conducting and 
insulating phases of interacting bosons in quenched disorder. In this effort to understand 
what is now known as the dirty-boson problem [5], most studies focused on the zero- 
temperature quantum phases, with a variety of approaches including Luttinger-liquid 
theory [6, 7], general scaling arguments [8, 9], Bogoliubov theory [10, 11, 12, 13], strong- 
coupling expansions [14], as well as numerical calculations with Monte-Carlo [15, 16, 17] 
and DMRG algorithms [18]. The picture that emerged revealed a rich interplay of 
bosonic statistics, disorder, repulsive interactions, and commensur ability effects in the 
presence of a lattice. The hallmark of this interplay is the restriction of the superfluid 
phase to a regime of moderate interactions and weak disorder, surrounded both at weak 
and strong interactions (or strong disorder) by a compressible gapless insulator called 
Bose glass [6, 9, 15, 16, 19]. To date, however, the characterization of the Bose-glass 
phase and the superfluid-insulator transition remains to a large extent an open problem. 
In spite of an important body of available results, this appears to be the case even for 
the one-dimensional geometry in view of recent publications [20, 21, 22, 23]. 
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The superfluid-insulator transition of disordered bosons attracted a renewed interest 
in the context of ultracold atoms [24], due to the high degree of control over disorder, 
interactions and conflning potentials achieved in these systems [25, 26, 27]. While the 
pioneering experiments with disordered bosons [28, 29, 30, 31, 32] aimed at observing 
Anderson localization in the noninteracting limit [33, 34, 35], more recent ones have 
provided flrst results towards a quantitative characterization of the phase diagram for 
nonvanishing interactions [36, 37, 38, 39, 40, 41, 42]. With this experimental activity, 
theoretical investigations also turned to the weakly interacting regime, which hitherto 
had been only poorly characterized. The scaling of the superfluid-insulator phase 
boundary as a function of the strength of disorder and interactions was established at 
the mean-fleld level [43, 44, 45, 46, 47], and shown to depend in an essential way on the 
microscopic disorder correlations. For the ID geometry, the fragmentation mechanisms 
driving the transition were analyzed by means of Bogoliubov theory [48] , while universal 
features of the transition and many-body corrections at intermediate disorder strengths 
were worked out with real-space renormalization group (RG) techniques [20, 49]. In the 
latter approach, making contact with experiments is a challenging task [49], and the 
method itself is not generalized to higher dimensions in a straightforward way [50]. 

To date, the details of the superfluid to Bose-glass transition in dimension d > 1 
are not well known, and the mechanisms driving the transition are expected to be more 
complex than in ID. The notion of weak links and fragmentation, for instance, involves 
connectivity in higher dimensions i.e. percolation [51, 50]. Moreover, while providing 
a natural step towards higher dimensions, the 2D case is particularly interesting in 
several respect. First, it stands for the marginal dimension of Anderson localization 
at the single-particle level, for the orthogonal and unitary Wigner-Dyson universality 
classes [52, 53]. Naively, one would therefore expect this geometry to be very sensitive 
to the introduction of interactions. Second, the clean weakly interacting system forms 
a true condensate at T = 0, and an algebraic superfluid for < T < Tbkt, where 
^EKT is the temperature of the Berezinksii-Kosterlitz-Thouless transition [54]. An 
outstanding question in this respect is how the zero-temperature superfluid-Bose glass 
transition connects to the clean BKT transition as disorder and temperature are varied. 
Experiments are ongoing in this regime to characterize the properties of the disordered 
Bose gas in 2D [55, 56]. 

The features of the 2D superfluid to Bose glass transition beyond mean fleld have 
been addressed recently. In Ref. [57], the Lifshits-tail physics associated with the deep 
insulating regime was analyzed by means of a multi-orbital Hartree-Fock method based 
on a set of low- lying single-particle states [57]. A real-space RG approach was devised 
for the 2D dirty-boson problem in Ref. [50] , and applied to the particle-hole-symmetric 
case, where the insulating phase is an incompressible Mott glass instead of a Bose 
glass. The analysis also emphasized the possible limitations of the strong-disorder RG 
in the 2D study. In Ref. [58] the T = phase diagram was studied by means of 
a weak-disorder expansion of Bogoliubov theory, valid far away from the transition, 
and quantum Monte-Carlo calculations. Upon flnite-size scaling of superfluid fractions 
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obtained numerically, the authors concluded in favor of a smooth crossover between the 
superfluid and insulating phases. However, the ad-hoc scaling law used in the analysis 
may deserve a careful examination, as the system sizes used in the numerics cover a 
relatively modest range. Hence, the development of a method reaching beyond mean 
field and affording large system sizes appears highly desirable in order to fill the existing 
gaps in the understanding of the two-dimensional case. 

In this contribution, we present a numerical scheme for the efficient computation 
of the one-body density matrix of a weakly interacting Bose gas in the framework 
of Bogoliubov theory, appropriately extended to account correctly for diverging phase 
ffuctuations in low dimensions [59]. The asymptotic behavior of the one-body density 
matrix determines the superfiuid or insulating behavior of the disordered Bose gas, and 
thereby allows a discussion of the phase diagram. Our scheme accomodates arbitrary 
disorder strengths, it has no intrinsic limitation in dimensionality, and it admits a 
straightforward extension to nonvanishing temperatures within the range of validity of 
Bogoliubov theory. The underlying density-phase representation allows for an accurate 
description of condensate, quasicondensate and insulating phases in the limit of large 
densities for any fixed interaction energy. The key feature of our approach is that it is 
based on an iterative scheme called the kernel polynomial method (KPM) [60], which 
allows the computation of correlation functions in large systems. 

The article is organized as follows. Section 2 provides a reminder on Bogoliubov 
theory in the density-phase formulation and on the form of the one-body density matrix 
within that framework. The KPM scheme for the computation of the one-body density 
matrix is detailed in section 3. In section 4, we validate our scheme by applying it to 
the case of disordered bosons at T = in ID and 2D, and by comparing our results with 
existing literature. In the ID geometry, we analyze the destruction of quasi-long-range 
order by disorder, and recover the phase diagram obtained through a direct approach 
in Refs. [48, 46]. In the 2D case, we compute the condensate depletion induced by both 
interactions and disorder, and compare our findings with those of Ref. [61]. In section 5, 
we conclude and discuss extensions of the present work. 



2. One-body density matrix of weakly interacting Bose gases 
2.1. Long-range and quasi-long-range order 

We consider a dilute gas of Bose particles described by the many-body Hamiltonian 

H = Jdr if\r)Hoif{r) + H\r)if\r)if{r)if{r) , (1) 

where ^^(r) is the bosonic field operator, 5^ > is the coupling constant parametrizing 
a repulsive contact interaction, and ^0 = + ^{^) the single-particle 

Hamiltonian. In the following, the external potential is a homogeneous random 
potential with zero mean, root-mean-square amplitude A, and a correlation length rj 
defined as the spatial width of the two-point correlation function V{r)V{r^). Here and 
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in the following, the bar denotes a statistical average over the disorder configurations, 
while the brackets (.) indicate quantum-mechanical expectation values. 

In the weakly interacting regime and close to the ground state, the properties of 
the dilute Bose gas are well described by standard Bogoliubov theory in 3D [62]. In this 
standard approach, the field operator ^!{r) is split into a classical component ^o('^) 
representing a condensate wave function with well-defined phase and a field 5^(r) 
describing quantum fiuct nations. An effective Hamiltonian is then derived from the 
expansion of H to second order in (^^. In dimensions < 2, however, the Mermin- 
Wagner-Hohenberg theorem [63, 64] rules out the presence of a condensate {i.e. long- 
range order [65]) at any temperature T > 0, as well as T = in ID. At sufficienty low 
temperatures in 2D and T = in ID, the Bose gas nevertheless forms a quasicondensate 
with a power-law decay of the one-body density matrix 

G(r,rO = {i!\r)i!{r')) (2) 

at large distances (||r' — r\\ +oc), and exhibits superfiuidity [21, 54]. While the 
absence of a true condensate in those cases precludes a straightforward application 
of standard Bogoliubov theory, the reformulation of the latter in a density-phase 
representation [66, 67] leads to a theory that is free of divergencies, and which captures 
the algebraic decay of correlation functions in quasicondensates [59, 68]. In this 
formulation the field operator writes ^(r) = e^^('^) ^p(r), where 6{r) is a phase operator 
that is safely defined in the high-density limit, and p(r) is the density operator. 
The latter is split into a classical component po(^) corresponding to the mean-field 
condensate (or quasicondensate) density and a fiuctuation Sp{r). At low temperatures, 
such an approach is accurate in the small-^' (large-po) limit for any constant U = gpo. 
This also holds in a disordered system, where the density po{r) may locally assume small 
values, but the regime of validity is recovered for sufficiently large at given U = g~p^. 

While the relation between superffuidity and condensation or quasicondensation 
is rather subtle, the presence of a compressible superffuid appears both necessary and 
sufficient for the existence of a condensate or quasicondensate [26]. Thus, the long- 
distance behavior of the one-body density matrix (2) allows a distinction between 
superffuid and insulating phases. In the clean (interacting) ID system, G(r, r') decays 
algebraically at T = 0. In the clean 2D system at T = 0, the one-body density 
matrix exhibits a plateau at long distances, which characterizes a true condensate. For 
< T < Tbkt, the 2D system is also an algebraic superffuid with a power-law decay 
of correlations. In 3D, ffnally, the system is a true Bose-Einstein condensate below the 
critical temperature Tbec- All these phases can be distinguished from the normal phases 
found at higher temperatures, which exhibit an exponential decay of G(r, r'). Similarly, 
the complete suppression of superffuidity by disorder at the phase transition to the Bose 
glass phase coincides with the destruction of long-range order or quasi-long-range order, 
i.e. with the emergence of an exponential decay of the one-body density matrix [46]. 
Since the ground state of the interacting Bose gas is globally extended and G(r,r') 
depends on the disorder conffguration, the disordered phases can be characterized by 
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the long-distance behavior of the statistical average 



mr-r'\\) = g,{ryi (3) 
where gi is the reduced one-body density matrix 

g (r r') - ^^^''^^^ (4) 
^p(r)p(rO' 



with p(r) the total gas density 

As anticipated above, G{r^r') and gi{r^r') can be calculated in a density-phase 
formulation of Bogoliubov theory [59]. Sections 2.2 and 2.3 provide a reminder on both 
the number-conserving and nonconserving variations of such a theory. These results 
are used in the following sections for the computation of the one-body matrix with the 
kernel polynomial method. 

2.2. Bogoliubov theory in number- conserving and nonconserving approaches 
In the ground state, the density po obeys the Gross-Pitaevskii equation (GPE) 



[i^o + gpo{r)] ^po{r) = fi^po{r), (5) 

where /x corresponds to the chemical potential in a grand-canonical description of the 
system. The quantum fluctuations and elementary (quasiparticle) excitations of the 
Bose gas are described by the fleld B{r) = 5p{r) / {2yJ~p^{r)) + i^ po{r)9{r)^ which 
obeys the equation of motion [59, 





= >Cgp At ■ (6) 



Here £gp is the standard Bogoliubov operator 

jC =( Hgp + gpoir) - ^ gpoir) . 
\ -gpoir) -Hgp - gpoir) + p 

with Hgp = Hq + gpoir). The field B admits the expansion [59, 70] 

Bir) = V L(r)5, + v*{r)b^] - i^oMr)Qsb + ^^6, (8) 



where bj is a bosonic quasiparticle operator, and the two last terms are explained below. 
The mode functions Uj{r) = {r\uj) and Vj{r) = {r\vj) are given by the solutions of the 
usual Bogoliubov-de Gennes equations (BdGEs) 

^-('i:;0^"'('i::0 

The operator jCqp is non-Hermitian, but its eigenvalues are real in the ground state of 
Hamiltonian (1). As /^gp* = ^gf^ the components Uj{r) = {r\uj) and Vj{r) = {r\vj) 
can be chosen as real- valued functions. We nevertheless consider the more general case 
of complex uj and Vj. For each eigenvector \vj))^ with eigenvalue Ej > 0, the 
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operator /^gp also has an eigenvector (|'^*), with eigenvalue —Ej < 0. The adjoint 

vectors of these positive and negative eigenvectors are —\vj))^ and (— 1^^*))^, 

respectively [71]. The biorthogonality of both positive and negative solutions of the 
BdGEs is thus expressed by the well-known relation 



(u,\uf) - (vM') = J drlu-(r)uAr) - v',(r)vAr)] = (10) 

The operator Cqf also has an eigenvector pertaining to the eigenvalue = 0, 
namely (|(/)o), — |(/>o))^, where (/)o{r) = ^y pQ{r) / Nq is the normalized ground state, with 
A^o = J drpo{r). As the set of eigenvectors of Cqf is not complete, an eigenvector 
(l^a), |0a))^ of with eigenvalue zero, such that (/>a('^) ^ K and ((/>o|0a) = 1/2^ is 
added to the set to obtain the closure relation [70] 

+ E ( ) (("il -("jD + ( III ) (-("ll' (".'D- (11) 

The Psb and Qsb terms in Eq. (8) account for fluctuations in the particle number 
and are responsible for the phase diffusion of condensates in symmetry-breaking 
approaches [72]. These terms do not arise in number-conserving approaches [73, 70, 74], 
which retain only fluctuations that are orthogonal to the ground state (/>o(t). The 
field A(r) that describes the orthogonal fluctuations obeys an equation similar to Eq. (6): 

(12) 

where Cqp has been replaced by [70] 

r=( Hgp + gQpo{r)Q - 1^ gQpo{r)Q \ , . 

\ -9Qpo{r)Q -HGP-9Qpo{r)Q + P )' ^ ' 

and (5 = 1 — |(/>o)(0o| projects orthogonally to the ground state. Equation (8) is replaced 
by the modal expansion 

Hr) = Y,uf{r)bj + vf{r)b], (14) 
where uj-{r) and Vj-{r) are solution of the modifled BdGEs 
\ui) \ _ / \u-^ 




J 



J E, y j^lj J iE, > 0). (15) 

The operator C has the same spectrum as jCqf^ and its positive and negative families of 
eigenvectors are simply obtained through the projections \uj-) = Q\uj) and \vj-) = Q\vj)^ 
which leave the biorthogonality relations (10) unaffected. Unlike jCqf^ however, the 
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operator C is diagonalizable. The zero eingenspace is spanned by the vectors (|(/>o), 0)^ 
and (0, |0o))^5 so that the resolution of identity writes [70] 

+ E ( ill ) -i^tn + ( I'lf-I ) (IS) 

As we shaU see below, Eqs. (5), (13) and (16) [or, equivalently, Eqs. (7) and (11)] are all 
that is needed for the efficient calculation of the one-body density matrix of disordered 
Bose gases in the weakly interacting regime. 



2.3. One-body density matrix in the density-phase representation 



The expression of the one-body density matrix G{r^r^) = (^^^(r)^(r')) was derived in 
the density-phase formalism in Ref. [59]. At zero temperature, it can be cast into the 
form [48, 75] 



G(r,r') = Jp{r)p{r') exp 



2 ^ 

j 



-(r) 



VMr) \/po{r') 



(17) 



where j enumerates the Bogoliubov modes with Ej > Q. This expression is valid in 
the limit of small density and phase fluctuations, which is realized at large average 
density for any given interaction energy U = g'po. In this limit, one has p(r) Po('^)5 
which in the presence of a true condensate amounts to a small condensate depletion (see 
section 4.2). In this limit also, the argument of the above exponential depends directly 
on the solutions of the GPE (5) and the BdGEs (15), without higher-order fluctuation 
corrections. 

Remarkably, the expression (17) accurately describes weakly interacting Bose gases 
in any dimension. In particular, it is not plagued by divergences in low dimensions, 
and captures the power-law decay of the one-body density matrix of quasicondensates 
of one-dimensional Bose gases at T = [59, 68, 76]. This expression was also used in 
Ref. [48] in conjunction with a numerical diagonalization of the Bogoliubov operator 
to analyze the destruction of quasi-long-range order by disorder in the ID geometry. 
While Eq. (17) involves all Bogoliubov modes, the sum is typically dominated by the 
modes of the low-energy phonon regime. However, even with a restriction to low-energy 
modes, the calculation of the one-body density matrix G(r,r') through complete or 
partial diagonalization of the Bogoliubov operator becomes prohibitive for large system 
sizes. In the following section, we present an alternative scheme, based on the kernel 
polynomial method [60], which circumvents the solution of the Bogoliubov eigenvalue 
problem and constitutes the main result of the present work. 
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3. Kernel polynomial scheme for the one-body density matrix 

The kernel polynomial method (KPM) [77, 78, 60] is an iterative numerical scheme for 
the computation of correlation functions. The KPM bypasses the spectral decomposition 
of the operators involved in those correlation functions, which may be numerically 
intractable. The KPM technique and some applications have been recently reviewed in 
Ref. [60]. In section 3.1 we introduce the elementary aspects of KPM that are relevant 
to the present study. In section 3.2, we show how a KPM scheme can be devised to 
compute the one-body density matrix G(r, r') on the basis of Eq. (17). 



3.1. Basics of the kernel polynomial method 

The core ingredient of KPM is the expansion of correlation functions on Chebyshev 
polynomials of the first kind. The latter are orthogonal polynomials on the interval 
/ = [—1,1], defined by the recurrence relation T^+i(x) = 2xT^(x) — T^_i(x) with 
To(x) = 1 and Ti{x) = x. Consider a Hermitian operator X with a discrete or continuous 
spectrum contained in /, and the correlation function 

/(|a),|6),x) = {a\6{X -x)\b), xeL (18) 

The latter is a formal writing for 

/(|a), |6),x) = ^S{xj - x){a\xj^a){xj,a\b), (19) 

where is an orthonormal eigenbasis of X, which provides the spectral 

decomposition X = ^^^^ with Xj G /, and the resolution of identity 

l = 5]ki,a)(2;„a|. (20) 

The above correlation function has the expansion 



f{\a),\b),x) ^ 



(21) 



M\a),\b)) + 2j2l^n{\a),\b))Tn{x) 

n=l 

where the coefficient /x^(|a), |6)), called Chebyshev moment of order n, is defined as 

f^n{\a),\b)) = J J{\a),\b),x)T^{x)dx. (22) 

Owing to Eqs. (19) and (20), the moments of /(|a), |6),x) are simply given by 

Mn(|a),|6)) = J2Tn{xj){a\xj^a){xj^a\b) = {a\Tn{X)Y^\xj^a){xj^a\b), (23) 

which boils down to the matrix element of a polynomial of X: 

/i„(|a),|6)) = (a|r„(X)|5). (24) 

Instead of calculating T^(X) for each new n index, which is computationally costly, one 
takes advantage of the recurrence relation between Chebyshev polynomials and keeps 
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track of two vectors, = T^(X)|6) and = Tn-i{X)\h)^ with the initiahzation 

|6o) = 1^) and = X\b). Then, a single apphcation of X (matrix- vector multiphcation) 
yields the new vector = 2X|6^) — l^n-i), along with the next Chebyshev moment 

/i^+i(|a), \h)) = {a\bn+i). In practice, the expansion (21) is truncated at some finite 
order A^, and /(|a), is approximated by 

1 



/(^)(|a),|6),x) 



N-l 



^rVo(la), \b)) + 2 V eVn(la), \b))Ux) 



n=l 



(25) 



where the g^"^ factors are introduced to damp the Gibbs oscillations caused by 
the truncation [60]. These factors amount to the action of a convolution kernel 
on /(|a), x), whence the name of KPM. Finally, the functional dependence of 
/^^^(|a), 1 6), x) on X is usually efficiently computed for a set of points x/^ G / by using a 
discrete cosine transform [60]. 

In summary, the KPM offers a simple iterative scheme for the calculation of 
correlation functions akin to expression (18), and avoids the numerical complexity 
associated with the spectral representation (19). Let us now examine how this method 
can be applied to the Bogoliubov operator for the computation of G(r, r'). 

3.2. Calculation of the one-body density matrix 

Expression (17) for the one-body density matrix involves the eigenmodes of the 
Bogoliubov operator the spectrum of which lies on the real line. As all subsequent 
numerical calculations are carried out with a finite-difference scheme and finite-size 
systems, the spectrum of C has a compact support [— £'max, ^max], where E'max depends 
on the tight-binding bandwith of the finite-difference scheme, the strength of interactions 
and the disorder configuration V{r). In the calculations presented in section 4, the upper 
bound E'max is obtained by solving the sparse Bogoliubov eigenvalue problem (9) for the 
largest eigenvalue with a Lanczos method. Taking a slightly larger E'max to ensure good 
KPM convergence at the spectrum boundaries [60], the spectrum is mapped to / by the 
rescaling C C/E^^y,. 

To exhibit correlation functions akin to Eqs. (18) and (19), we cast expression (17) 
into the following form: 



G(r,r') = Wp(r)p(r')exp 



1 

- F{r,r\e)de 



(26) 



with 



and 



F(r, r', e) = /(r, r, e) - /(r, r', e) - /(r', r, e) + /(r', r', e) (27) 



/(r, r , e) = - 2^ 5{ek - e) — 



(28) 
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where {w^^^) |t^^2) ^he second components of the eigenvector {{wk^i)-, |'^/c,2))^ of 
and its adjoint vector {\w^i)^ \^k2)Y ^ respectively. In Eq. (28), the index k runs over 
the positive [E^ > 0), nuU (Ej. = 0) and negative {Ej. < 0) famihes of eigenvectors, and 
Ck = £'/e/£'max- Bccausc of the integration boundaries in Eq. (26), the term f{r^r^^e) 
contributes to the exponent by 



r VPoir)po{r'y 



(29) 



where j enumerates the modes with Ej > 0. The —1/{2Nq) term, which stems from 
the zero eigenvector (0, l^o))^, cancels out of the / sum in Eqs. (26) and (27), so that 
Eq. (26) is indeed equivalent to Eq. (17). 

The modes with Ek < are included in sum (28) to use the resolution of identity 
[see Eq. (16)] 

1 = E(Ki)' K2)f{{Kil K2I). (30) 
k 

Indeed, rewriting Eq. (28) as 

(0, {r\){\wi,), H,)n{wi,i {wy)io, y)r 



we find that the Chebyshev moments of /(r, r', e) are 



(31) 



and those of F{r^r\ e) follow as 

Mn{r, r') = iin{r, r) - jinir, r') - /i^(r', r) + /i^(r', r'). (33) 

Finally, there is no need for a Chebyshev inversion by discrete cosine transform, as the 
expansion (21) can be integrated analytically on [0, 1], and we obtain 

^ ' ' ^ 2 TT^ 2p+l ^ ^ 

Note that the contributions of all even moments except Mo(r, r') are integrated out on 
[0, 1]. The moments /x^(r,r') are calculated iteratively following the recurrence scheme 
outlined in section 3.1, for the four (r,r') pairs in M^(r,r'). This requires only two 
Chebyshev sequences as, for instance, r^(>C/£'max)(0, be used to compute 

/i^(r,r') and /x^(r',r'). When the Chebyshev iterations are truncated at order A^, the 
reduced one-body density matrix is approximated by 



exp 



^p(r)p(r') 



(N) Lt-1J (_l)pJN) 



(35) 



where the g^"^ are convolution kernel factors (see section 3.1). We found the standard 
Jackson kernel [60] to be suitable in this scheme. 
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The Chebyshev iterations based on Eq. (32) involve projections orthogonally to the 
ground state by means oi Q = 1 — |(/)o)((/>o|. Interestingly, the Bogoliubov operator C 
may be replaced by /^qp in the iterations, so that projections are not necessary. This 
provides a further simplification of the KPM calculation of the one-body density matrix. 
Indeed, upon expansion with vj-{r) = Vj{r) — ((/)o|'^j)(/>o(T), one easily sees that Eq. (17) 
also writes as 

2" 



G(r,r') = Wp(r)p(r')exp 



Vpoir) \/po{r') 



(36) 



The function FGp(r, r',6) is defined as the analog of F{r^r^^e) in Eqs. (26) and (27), 
with four terms of the form 

{r\Wk'^2){'Wk',2W) 



fGp{r,r',e) = -^5(6^/ - e)- 



(37) 



fc/ ^^0(^)^0(0 

where k' enumerates eigenmodes of Cqf with Ej.^ ^ 0. Owing to the closure relation (11), 
the Chebyshev moments of /gp(^, e) read 



where 



yP^yr) V^max/ VPOI'T ) 



0o(r-') 



E 

>^GP 

J-^ma.-K- 



Given that £gp(|<^o), -|0o))^ = (0,0)^ and £Gp(|0a), |0a))^ = «(|<^o), 
constant a, as detailed in Ref. [70], we find 

R^iry) ={-iy 



R^iry) ={-iy 



Mr) 



y/Nopo{r) 



fi?^(r,r') = (-lF«<|^. 

^Vo^max 



(39) 

(40) 
bo))'^ with 

(41) 

(42) 
(43) 

(44) 



These Rn\r^r^) terms cancel out of the sum 

M^^{r,r') = ^,^^iry - ^,^^{ry) - f,^^y,r) + ^,^^y,r'). (45) 

Hence, the comparison of Eqs. (32) and (38) shows that Cqp can be used instead of C 
in the Chebyshev iterations underlying Eq. (35). 

The expressions (32), (33) and (35) are the main results of this section. They 
provide an iterative scheme for the calculation of the reduced one-body density matrix of 
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a weakly interacting Bose gas, once the solution of the GPE (5) is given. In the following 
section, this scheme is applied in various geometries, with and without disorder. 

Remarkably, the above approach can be extended in a straightforward way to 
nonzero temperatures within the framework of Bogoliubov theory. In addition to 
the Vj{r) terms, the thermal G{r^r') also contains contributions from the Bogoliubov 
components Uj{r)^ which can be calculated through additional KPM iterations. As the 
Uj and Vj terms are weighted by Bose-Einstein occupation factors, the integration in 
Eq. (34) no longer has a simple analytical solution. However, this analytical step may 
be replaced by a discrete cosine transform at low computational expense. 

Our results also show that the operator involved in the correlation function and 
the KPM iterations need not be Hermitian. Some illustrations of this fact can be 
found in the literature, with special cases such as the computation of retarded Green's 
functions [60] or the solution of generalized Hermitian eigenvalue problems [79]. The 
Bogoliubov operator C and >Cgp provide two other interesting examples in that context. 
First, C is diagonalizable, albeit non-Hermitian, and its eigenvalues are real. As a 
consequence, the eigenvectors and their adjoint s form a complete biorthogonal set that 
can be used for the closure relation, and there is no need for a twofold KPM iteration to 
retrieve the spectral information on a compact set of the complex plane. Second, Cq^ is 
not diagonalizable, and yet its generalized eigenvectors (and their adjoints) can be used 
for the closure relation. In the derivation of Eqs. (41) to (44), we took advantage of 
the fact that (|(/)o), — |(/>o))^ and (|(/>a), \4^a)Y generalized eigenvectors of low order, 
so that the action of T^(i3Gp) can be evaluated easily. 



4. Application to disordered bosons 

We now employ the KPM scheme introduced above to analyze the asymptotic behavior 
of the one-body density matrix of disordered bosons in ID and 2D. While the approach 
of the previous sections is general, we consider here a Gaussian random potential V 
with Gaussian correlation function 



V{r)V{r') = A^e-^'--'-^ ^ (46) 

The numerical procedure for the computation of gi{r^ r') is the following. The continuum 
problem (1) is represented on a finite volume in a finite-diflFerence scheme with lattice 
spacing I = L/n^. To emulate the continuum limit, the hopping term t = fi^/{2mP) 
is chosen to be much larger than all the other energy scales of the problem. In all the 
calculations presented here, the correlation length of the disorder is taken to be = 4^, 
which is sufficient for our purposes. The correlation length 77 and the associated energy 

^ 2mrf ^ 
are used as reference scales throughout this section, even in the absence of disorder 
(A = 0). For each configuration V^(r), the ground-state solution Po{r) of the 
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GPE (5) and the corresponding chemical potential /x are computed through imaginary- 
time propagation with a standard Crank-Nicholson scheme, at fixed particle number 
No = J drpo{r). We denote by 

(48) 

the average mean-field interaction energy. Periodic boundary conditions are imposed on 
the GPE and on the BdGEs (9). The Bogoliubov operator Cqf is set up on the basis 
of V{r) and po{r). Then, gi{r^r') may be calculated with the KPM iteration detailed 
in the last section for any (r, r') pair. 

4^1. Superfluid to Bose glass transition in one dimension 

In Refs. [48, 46], the destruction of quasi-long-range order was used as a signature of the 
superfiuid to Bose glass transition at T = in ID, and this criterion was used to draw 
the quantum phase diagram on the basis of the asymptotic behavior of the (averaged) 
reduced one-body density matrix ^(|r — Here, we use the ID setting as a testbed 
for the KPM approach described above. 

In the absence of disorder, gi is expected to decay at large distances with a power 
law given by [80, 59]: 

g^{r,r')^(- ^ , |r-/|>e, (49) 



^4|r — 

where C = 0.57721 ... is Euler's constant, the density p can be approximated by po 
within our Bogoliubov approach, and ^ is the healing length defined here as 

with U = gpo. Figure 1 shows the result of a KPM calculation for a system of 
length L = 2^^7y, and the excellent agreement obtained with the power law (49) for 
\r — r^\ > i.e. in its regime of validity. The regrowth observed at large \r — r^\ is 
due to the periodic boundary conditions, and does not aflFect significantly the data for 
\r — r^\ < L/4:. In all the subsequent analyses, we retain only this range for determining 
the asymptotic behavior of gi. Note the large system size achieved in this computation. 
In this homogeneous case, a single KPM iteration suffices for the computation of gi{r^ r'). 
The number of moments required to resolve all individual Bogoliubov modes in the 
phonon regime and achieve a good convergence of the KPM result grows linearly with 
the system size = L/£. The required storage space is also of a few n^. This has to be 
compared to the storage space and computation time required for a full diagonalization 
of >Cgp, which scale as nj and respectively [60]. 

For A > 0, the reduced one-body density matrix gi{r^r') depends on the disorder 
configuration, and is no longer translation invariant. Fig. 2 shows the behavior of 5^1(0, r) 
for a single disorder configuration, and the results obtained with the KPM scheme for 
various moment numbers [see Eq. (35)]. When is sufficiently high, the KPM results 
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Figure 1. Reduced one-body density matrix gi{r^r') in the absence of disorder for 
a ID system of length L = 2^^?7, with interaction strength U = 0.10£^c (^•6- ^ — 14?7). 
The blue solid line is the result of the KPM calculation. The red dashed line is the 
asymptotic power law given by Eq. (49) with p pQ. 

coincide with those obtained from complete diagonahzation, and faithfuUy reproduce the 
spatial details of gi{r^r'). As a general trend, we also observe that the KPM estimates 
of gi{r^r') converge from above with increasing A^. This can be attributed to the fact 
that low N values imply poor spectral resolution, and hence are not able to resolve the 
small energy scales associated with the long-distance decay of gi , such as for example the 
vanishing energy separation that can arise when Bogoliubov modes are strongly localized 
in different spatial regions. While the number of moments required for convergence 
scales linearly with the system size in the clean case, we also observed that this scaling 
is slightly faster than linear in the disordered case. 

After averaging over disorder, the one-body density matrix exhibits either a power- 
law or an exponential decay a large distances, depending on the strength of interactions 
and disorder. For long-enough systems a similar behavior may already be observed 
qualitatively with a single disorder configuration in the spatial average 

1 

Qiir) = -J gi{r',r' + r)dr'. (51) 

Fig. 3 shows spatial averages computed with the KPM scheme for two sets of parameters 
in a system of length L = 512r]. We display here the range r ^ L/A^ and the crossover 
to the long-distance behavior is visible on both panels. The linear behavior found for 
U = lAAEc on the double logarithmic scale indicates a power-law decay of the averaged 
g'l, corresponding to the superfluid phase. For U = 0.48£'c, on the other hand, we find 
an exponential decay indicating a Bose glass. For the same system size, the convergence 
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Figure 2. Reduced one-body density matrix ^i(0,r) for a single disorder 

configuration, in a ID system of length L = 128r]. The interaction strength is 
U = gNo/L = O.SEc, and the disorder strength is A = 0.7 Ec. The black dashed 
line shows the result obtained from the complete diagonalization of the Bogoliubov 
operator jCqp- The other curves are the KPM results obtained for various numbers of 
moments. 
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Figure 3. Spatial average (r) for a single disorder configuration with amplitude 
A = O.SEc and two interaction energies /7, in a ID system of length L = 512r]. The case 
U = 1.12Ec (left panel) exhibits a power-law decay of , while the case U = 0A8Ec 
(right panel) exhibits an exponential decay for r > 6O77. 



of the KPM result in the insulating phase requires a higher number of moments than 
in the superfluid phase. This may be attributed to the increase of the Bogoliubov 
density of states near zero energy and to a reduced level repulsion. Indeed, as the 
disorder increases, the system turns progressively into a collection of superfluid puddles 
separated by high potential barriers. 
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Figure 4. Ground-state phase diagram of weakly interacting disordered bosons in ID. 
The phases are characterized by the power-law decay (blue circles; superfluid) or the 
exponential decay (green squares; Bose glass) of the averaged one-body density matrix 
gi{r). The latter was obtained through KPM iterations with system sizes varying 
between 64?] and 204877. The black triangles lie on the estimated phase boundary. 
Linear fits to these data points on the double logarithmic scale yield the slope 0.75±0.04 
in the white- noise regime U <^ Ec, and 0.94±0.03 in the Thomas- Fermi regime U ^ Ec. 

The quantum phase diagram of the weakly interacting regime can be drawn by 
varying A and [/, and determining for each parameter set whether the asymptotic part 
of the disorder-averaged gi behaves as a power law or an exponential This procedure 
is put on a systematic footing by fitting the long-distance part of ^(r) with both 
power laws and exponentials, and monitoring the fit quality. Fig. 4 shows the phase 
diagram obtained with the KPM technique. The blue circles and green squares have 
been identified as belonging to the superfluid and Bose glass phases, respectively. The 
black triangles cannot be attributed to one phase or the other with the accuracy of the 
data, and are assumed to lie on the phase boundary. The red solid lines are fits to the 
black triangles in the regimes U % and U ^ E^. In the white-noise regime U <^ E^ 
{i.e. ^ ^ ry), the phase boundary is expected to scale as A/Ec ^ (U/Ec)^/^^ while in the 
Thomas- Fermi regime U ^ E^ {i.e. C ^ v) the critical A is expected to grow linearly 
with U [48, 45, 46, 81]. The fits of Fig. 4 are in good agreement with these scaling laws. 
Our findings thus reproduce the results of Refs. [48, 45] without the need for partial or 
complete diagonalization. This validates our approach. 
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4^2. Condensate depletion in two dimensions 

In the absence of disorder, weakly interacting bosons form a true condensate at T = 
in 2D. In this case, the reduced one-body density matrix gi{r^ r^) tends to a constant for 
||r — r'll ^ oo, equal to the condensate fraction po/p- To leading order in the strength 
of interactions, the condensate fraction is given by [82, 54, 59, 61] 

where g is the 2D coupling constant. While scattering theory shows that this constant 
actually depends logarithmically on the 2D density p and the 3D scattering length a [54] , 
we take it as a given parameter of Hamiltonian (1) in 2D, even for the strongly 
inhomogeneous case. According to Eq. (52), interactions reduce the condensate fraction. 
In the presence of disorder, the condensate fraction is expected to be further reduced 
(but nonzero) in the superfluid regime, and completely suppressed in the Bose glass 
phase, due to the (exponential) decay of the one-body density matrix. 

The KPM algorithm introduced above is expected to reduce significantly the 
computational cost of a study of the superfiuid-insulator transition on the basis of 
the one-body density matrix. A fuUy-fiedged analysis of this transition nevertheless lies 
beyond the scope of the present paper, and is left for future work. Here, we restrict 
ourselves to the superfiuid regime and consider the calculation of the disorder-induced 
condensate depletion as a benchmark for the KPM technique. This depletion has been 
calculated analytically in Ref. [61] for the limit of weak interactions and weak disorder. 
To leading order in A/[/, the depletion 6p = p{r) — po{r) reads 



Sp 5p 



(0) 



(53) 



where 5p^^^ is the interaction-induced condensate depletion in a clean system, given by 
Eq. (52). The function h depends only on the ratio of the disorder correlation length rj 
and the healing length ^. Note that h differs from the similar function introduced in 
Ref. [61] by a trivial factor ^/2 due to a diflFerent definition of ^. 

The left panel of Fig. 5 shows the result of a KPM calculation of 5^1(0, r) in a clean 
2D system. Starting from the origin, the one-body density matrix drops and reaches 
a plateau beyond a few healing lengths. The regrowth of gi on the system edges is 
due to the periodic boundary conditions, as in the ID case. The condensate fraction 
is extracted from the value assumed at the center of the system. With this procedure, 
we studied the dependence of the condensate fraction on the interaction strength. The 
numerical results are plotted in Fig. 6, and stand in perfect agreement with Eq. (52). 

The disordered case was examined with a similar procedure. The right panel of 
Fig. 5 displays the values of 5'i(0,r) obtained for a single disorder configuration with 
A/[/ = 0.125. For this value the effect disorder is weak, and the one-body density matrix 
still exhibits a plateau, at roughly the same level than the clean case shown in the left 
panel. The statistical average and fiuctuations of the condensate fraction are evaluated 
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Figure 5. Reduced one-body density matrix gi (0, r) for a 2D system of size 64?] x 64/^, 
at T = 0, in the absence of disorder (left panel) and for a single disorder configuration 
with amplitude A = 4.0^c (right panel). The interaction strength is [/ = 32.0^c in 
both cases. 




Figure 6. Condensate fraction po/p versus interaction strength for a 2D system with 
a particle density of Nq/L'^ = 160r]~^. The abscissa parameter is g' = 2mg/h^. The 
blue circles are extracted from the values of the one-body density matrix in systems of 
size 64/^ X 64/^. The red line represents Eq. (52). 

by calculating 5'i(0,r), with r at the system center, for several disorder configurations. 
Figure 7 shows the average and the standard deviation of the fractional depletion 5p/ p 
[i.e. the complement of the condensate fraction) as a function of the disorder strength. 
For weak disorder, we indeed find a good agreement with the theoretical prediction (53), 
as shown in the inset. For A > 1.3£'c {i.e. A > 0.2[/), however, the numerical averages 
clearly lie below the theoretical curve. While a precise comparison at strong disorder 
would require a distinction between the local fluctuations of the depletion and those of 



Superfluid-insulator transition in weakly interacting disordered Bose gases 



19 




°'^0 5 10 15 20 25 30 35 

Figure 7. Fractional condensate depletion versus disorder strength in a system 
of size 64?] x 64/^, for an interaction energy U = 6AEc. The numerical data in 
blue represent the statistics of the reduced one-body matrix, evaluated by KPM at a 
distance ^/2L/2 (see text) for 100 to 200 disorder configurations. The dots correspond 
to average values, and the error bars indicate the root-mean-square fiuctuations around 
those averages. The red solid line represents the weak-disorder prediction (53) for the 
depletion, normalized by the average density. 

the fractional depletion, this suggests that higher orders in the weak-disorder expansion 
should reduce the condensate depletion. A detailed study of this effect lies beyond the 
scope of this work, as the analytical calculation of higher-order corrections becomes 
rather cumbersome. 

5. Summary and outlook 

We have developed an iterative scheme, based on the kernel polynomial method, for the 
efficient computation of the one-body density matrix of weakly-interacting Bose gases 
in the framework of Bogoliubov theory. Such a scheme is relevant for regimes of strong 
disorder, which cannot be tackled analytically. The scheme was applied to the case 
of disordered bosons at T = in one and two dimensions. In the one-dimensional 
case, we characterized the superfluid-insulator phase transition on the basis of the 
long-range behavior of the one-body density matrix, and successfully reproduced the 
results of Refs. [48, 46] with a low computational overhead. In the two-dimensional 
geometry, we analyzed the quantum depletion induced by interaction and disorder 
in the superfluid regime, and found a good agreement with results available for the 
weakly disordered regime. These case studies validate our approach and suggest that it 
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may be used to study the coherence properties of weakly interacting Bose systems for 
system sizes that remain hardly tractable with other numerical techniques. This feature 
is particularly interesting for investigations into the superfluid-insulator transition in 
higher dimensions. As outlined here, our approach is also easily extended to regimes of 
low but nonzero temperatures, which are relevant to ongoing experiments with ultracold 
atomic gases. 
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